Использовалась библиотека stats

Библиотеки

library(Rssa)
library(ggplot2)
library(mFilter)
data(USUnemployment)

Периодограммы

к 14.03

Эффект растекания синуса

Растекание — эффект, возникаемый, если длинна ряда не будет кратна периоду. Найдем длинну нашего ряда:

w <- 0.1  # Тут задаем частоту, которая соответствует периоду T = 10
x <- 1:100
y <- sin(2 * pi * w * x)
length(y)
## [1] 100

Получаем длинну ряда 100, построим периодограмму для ряда у которого длинна кратна периоду:

sp <- spec.pgram(y, spans = NULL, detrend = FALSE, log = "no", fast = FALSE, taper = 0)

Видим, что эффекта растекания на периодограмме нет, пик единственный и соответсвует заданной частоте. Теперь попробуем изменить длинну исходного ряда. Можно добавить несколько нулевых наблюдений в ряд, но лучше убрать уже имеющиеся, чтобы не искажать картину нулями:

yCut <- y[1:(length(y) - 3)] # Убираем три наблюдения
length(yCut)
## [1] 97

Теперь длинна ряда 97, что не кратно периоду T = 10, теперь посмотрим на его периодограмму:

sp <- spec.pgram(yCut, log = "no", detrend = FALSE, fast = FALSE, taper = 0)

Видим, что пик начал расползаться по соседним частотам.

Генерация белого шума

Шум — это случайная составляющая временного ряда. Если средний вклад всех частот одинаковый, то шум называется белым. Свойство периодограммы белого шума: значение периодограммы имеет экспоненциальное распределение с одним и тем же средним.

n <- 1000
wnoise <- rnorm(n, 0, 1)
spec.pgram(wnoise, detrend = FALSE, log = "no", fast = FALSE, pad = FALSE,taper = 0)

Так как для шума вообще нет детерминирующих частот, мы видим из периодограммы, что каждая частота входит примерно с одинаковой силой. С увеличением количества наблюдений будет увеличиваться значение каждой частоты на сетке:

nLarge <- 10000
wnoiseLarge <- rnorm(nLarge, 0, 1)
spec.pgram(wnoiseLarge, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE, taper = 0)

Генерация красного шума

Отличие красного шума от белого в том, что предыдущие значения влияют на следующие, т.е. имеют ненулевую корреляцию, в отличие от белого шума.

\(r_n = cor \cdot r_{n-1}+\sqrt{(1-cor^2)} \cdot \omega\)

w0 <- wnoise[1]
wnoise <- wnoise[2:n]
cor <- 0.8
rnoise = Reduce(function(prev_v, next_v) cor * prev_v + next_v * sqrt(1 - cor^2), wnoise, w0, accumulate = T)

spec.pgram(rnoise, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

Построение периодограммы для реальных данных

График безработицы в США

О данных: https://link.springer.com/chapter/10.1007/978-1-4612-5098-2_66 408 ежемесячных наблюдений безработицы в США среди мужчин Длину окна берем равную 12 в соотвествии с годовой периодичностью. Посмотрим на их график:

library(Rssa)
data(USUnemployment)

unempl.male <- USUnemployment[, "MALE"]

plot(unempl.male, type = "l")

Периодограмма безработицы в США

По графику можно сказать, что временной ряд достаточно волатильный, имеет сложный тренд, по периодограмме можно сделать выводы о частотах входящих в него:

spec.pgram(unempl.male, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

Видим, что наибольший вклад имеют низкие частоты, которые мы относим к шуму, также есть пики в точках w = 1/12 и w = 2/12, что говорит нам о наличии годовой и полугодовой сезонности.

Фильтры

к 14.03

Функция для фильтра скользящего среднего

movingAverage <- function(x, n=1, centered=TRUE) {
      
      if (centered) {
            before <- floor  ((n-1)/2)
            after  <- ceiling((n-1)/2)
      } else {
            before <- n-1
            after  <- 0
      }

      s     <- rep(0, length(x))
      count <- rep(0, length(x))
      
      new <- x
      count <- count + !is.na(new)
      new[is.na(new)] <- 0
      s <- s + new

      i <- 1
      while (i <= before) {
            new   <- c(rep(NA, i), x[1:(length(x)-i)])
            
            count <- count + !is.na(new)
            new[is.na(new)] <- 0
            s <- s + new
            
            i <- i+1
      }
      
      i <- 1
      while (i <= after) {
            new   <- c(x[(i+1):length(x)], rep(NA, i))
            
            count <- count + !is.na(new)
            new[is.na(new)] <- 0
            s <- s + new
            
            i <- i+1
      }
      s/count
}

АЧХ

afc <- function(filter, omega)
{
      k <- seq_along(filter) - 1
      h <- function(o) sum(rev(filter)*exp(-k*1i*o))
      abs(sapply(omega, h))
}

Применение фильтра скользящего среднего к модельным данным

Задаем модель

n = 100
x <- 1:n
period <- 10

modelTs <- (2*x + 3) +
 5*cos(2*pi*(1/period)*x) + 5*cos(2*pi*(2/period)*x) + 5*cos(2*pi*(3/period)*x) + 5*cos(2*pi*(4/period)*x) +
  + rnorm(n)

plot(modelTs, type = "l",  main = "Исходный временной ряд")

Исходная периодограмма

sp <- spec.pgram(modelTs, detrend = FALSE, log = "no", fast = FALSE, pad = FALSE, taper = 0)

Применяем фильтр

plot(modelTs, main = "Исходный ряд + Cкользящее среднее", type = "l")
legend(x = "topleft",                                # Position
       legend = c("timeseries", "MA(5)", "MA(10)"),  # Legend texts
       lty = c(1, 1, 1),                             # Line types
       col = c("black", "red", "blue"),              # Line colors
       lwd = 2)                                      # Line width

maHalf <- movingAverage(modelTs, 5, TRUE)
lines(maHalf, col = "red")

maFull <- movingAverage(modelTs, 10, TRUE)
lines(maFull, col = "blue")

Периодограмма MA(5)

sp <- spec.pgram(maHalf ,detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

sp <- spec.pgram(maHalf ,detrend = FALSE, log = "yes",fast = FALSE, pad = FALSE,taper = 0)

Периодограмма MA(10)

sp <- spec.pgram(maFull ,detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

sp <- spec.pgram(maFull ,detrend = FALSE, log = "yes",fast = FALSE, pad = FALSE,taper = 0)

АЧХ фильтра скользящего среднего

freq <- seq(0, pi, 0.001)
filt <- rep(1/period, period)
omega <- freq/2/pi

str <- c("0", "1/10", "2/10", "3/10", "4/10","5/10")
plot(afc(filt, freq) ~ omega, type = "l", xlab = "Frequency", ylab = "Frequency response", xaxt = "n")
title(main = "АЧХ фильтра скользящего среднего MA(10)")
axis(1, at = seq(0, 0.5, by = 1/10), las = 2, labels = sprintf("%s", str))

period = 5
filt <- rep(1/period, period)

plot(afc(filt, freq) ~ omega, type = "l", xlab = "Frequency", ylab = "Frequency response")
title(main = "АЧХ фильтра скользящего среднего MA(5)")

Применение фильтра к реальным данным

График безработицы в США Периодограмма безработицы в США

Для выделения тренда выставим длинну окна равную периоду ряда (T = 12). Наложим график скользящего среднего на исходный:

plot(unempl.male, main = "Исходный ряд + Cкользящее среднее", type = "l")
legend(x = "topleft",                                
       legend = c("timeseries", "MA(12)", "MA(24)"),
       lty = c(1, 1, 1, 1),                           
       col = c("black", "blue", "red"),    
       lwd = 2)                                     

ma12 <- movingAverage(unempl.male, 12, TRUE)
ma24 <- movingAverage(unempl.male, 24, TRUE)
lines(ma12, col = "blue")
lines(ma24, col = "red")

Теперь построим периодограмму скользящего среднего:

spec.pgram(ma12, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

SSA

к 14.03

Выделение тренда с помощью SSA

Применим SSA к реальным данным, ширину окна можно взять L = K/2 и кратную периоду (условие ассимптотической разделимости). Построим матрицу взвешенных корреляций, чтобы проанализировать компоненты:

Матрица взвешенных корреляций

tsSSA <- ssa(unempl.male, L = length(unempl.male) %/% 2)
tsSSACor <- wcor(tsSSA, groups = 1:30)
plot(tsSSACor)

Видим, что после 4 компоненты начинается зашумленность, при этом 15-я компонента достаточно хорошо отделяется от соседних.

Собственные вектора

plot(tsSSA, type="vectors", idx = 1:20)

По парам собственных векторов можно сказать о векторах периодичностей:

plot(tsSSA, type="paired", idx = 2:20)

Видим правильные фигуры в компонентах (5, 6) и (12, 13):

tsSSASeason <- reconstruct(tsSSA, groups = list(season = c(5, 6, 12, 13)))
spec.pgram(tsSSASeason$season, log = 'no', fast = FALSE, taper = 0, detrend = FALSE)

Видим, что выделились пики соответсвующие годовой и полугодовой периодичности, как на исходной периодограмме.

tsSSASignal <- reconstruct(tsSSA, groups = list(c(1:4, 7:10), c(5, 6, 12, 13)))

plot(tsSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))

FOSSA

Алгоритм улучшает сильную разделимость.

FOSSA: матрица корреляций

tsFOSSA <- fossa(tsSSA, nested.groups = c(1:13))
tsFOSSACor <- wcor(tsFOSSA, groups = 1:30)

plot(tsFOSSACor)

FOSSA: собственные вектора

plot(tsFOSSA, type = "vectors", idx=1:15)

FOSSA: итоговое разложение

tsFOSSASignal = reconstruct(tsFOSSA, groups <- list(c(5:13), 1:4))

plot(tsFOSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))

Корреляции SSA vs. FOSSA

plot(tsSSACor)

plot(tsFOSSACor)

Сезонности SSA vs. FOSSA

matplot(data.frame(tsSSASignal$F2, tsFOSSASignal$F2), type = 'l', col=c("red","green"), lty=c(1,1))

Тренд SSA vs. FOSSA

matplot(data.frame(c(unempl.male), tsSSASignal$F1, tsFOSSASignal$F1), type = 'l', col=c("black", "red","green"), lty=c(1,1))

Различные способы выделения тренда

к 21.03, 28.03

Тренд по POLY

trend.poly3 = lm(unempl.male ~ poly(1:408, degree = 3))
trend.poly5 = lm(unempl.male ~ poly(1:408, degree = 5))
trend.poly7 = lm(unempl.male ~ poly(1:408, degree = 7))


plot(c(unempl.male), main = "Исходный ряд + POLY", type = "l")
legend(x = "topleft",                                
       legend = c("timeseries", "POLY(3)", "POLY(5)", "POLY(7)"),
       lty = c(1, 1, 1, 1),                           
       col = c("black", "red", "blue", "green"),    
       lwd = 2)    

lines(trend.poly3$fitted.values, col="red", type = "l")
lines(trend.poly5$fitted.values, col="blue", type = "l")
lines(trend.poly7$fitted.values, col="green", type = "l")

Тренд по LOESS

LOESS метод основан на сглаживании с помощью построения локальных взвешенных линейных регрессий

index <- 1:408

fit.loess <- loess(unempl.male ~ index, span = 0.2, degree = 1)      
trend.loess <- predict(fit.loess)

plot(c(unempl.male), main = "Исходный ряд + LOESS", type = "l")
legend(x = "topleft",                                
       legend = c("timeseries", "LOESS"),
       lty = c(1, 1, 1, 1),                           
       col = c("black", "blue"),    
       lwd = 2)                                     
lines(trend.loess, col = "blue", type = "l")

Тренд по HP

chest = hpfilter(unempl.male, freq=60000, type="lambda")
trend.hp = chest$trend
plot(unempl.male, type="l")
lines(trend.hp, col="red")
legend("topleft", legend = c("Real TS","HP"), col = c("black","red"), lty=1)

Итоговый график со всеми способами выделения тренда

plot(unempl.male, type="l")
lines(trend.hp, col="red")
lines(ma24, col="blue")
lines(ts(trend.loess, start=c(1948, 1), frequency=12), col="orange")
lines(ts(trend.poly7$fitted.values, start=c(1948, 1), frequency=12), col="green")
legend("topleft", legend = c("Real TS","HP", "MA(24)", "LOESS", "POLY(5)"), col = c("black","red", "blue", "orange", "green"), lty=1)

Сравнение методов разложения

04.04

STL

STL имеет следующие параметры:

  • (inner) – число итераций внутреннего цикла (р).

  • (outer) – число итераций внешнего цикла. У нас в данных нет аутлаеров, поэтому outer = 0

  • (l.window) – сглаживающий параметр для low-pass фильтра. Возьмем равным 13. Рекомендуют брать ближайший нечетный к периоду.

  • (t.window) – сглаживающий параметр для тренда. Рекомендуют брать ближайшим нечетным к (1.5*period) / (1-(1.5/s.window)). Равен 23.

ns - (s.window) – сглаживающий параметр сезонности. Берем 13, так как годовая периодичность и параметр должен быть нечетным

tsSTL <- stl(unempl.male, s.window = 13, l.window = 13, outer = 0, inner = 1, t.window = 23)
plot(tsSTL)

Classical seasonal decomposition

tsCSD <- ssa(unempl.male, L = length(unempl.male) %/% 2, force.decompose =  FALSE, svd.method = "nutrlan")
fit.dec <- decompose(tsCSD, neig = 40)
tsCSDCor <- wcor(fit.dec, groups= 1:40)
plot(tsCSDCor)

plot(fit.dec, type = "vectors", idx = 1:20)

tsCSDSignal <- reconstruct(fit.dec, groups = list(trend = c(1:4, 7:10, 15), season = c(5, 6, 12, 13)))
plot(tsCSDSignal)

Итоговый график трендов разных разложений

plot(unempl.male, type="l")
lines(tsSSASignal$F1, col="red")
lines(tsFOSSASignal$F1, col="blue")
lines(tsCSDSignal$trend, col="orange")
lines(ts(data.frame(tsSTL$time.series)$trend, start=c(1948, 1), frequency = 12), col="green")
legend("topleft", legend = c("Real TS","SSA", "FOSSA", "CSD", "STL"), col = c("black","red", "blue", "orange", "green"), lty=1)

Примеры применения Toeplitz SSA и SSA with projections

SSA with projections

Projection SSA применяется для улучшения разделимости линейного тренда и сезонности. Используем смоделированные данные с линейным трендом (Применение фильтра скользящего среднего к модельным данным):

x <- 1:100
modelTs <- (2*x + 3) +
 + 5*cos(2*pi*(1/25)*x) + cos(2*pi*(1/5)*x) + rnorm(50)

plot(modelTs, type="l")

# Применим SSA с двойным центрированием для лучшего отделения линейного тренда
modelTsSSA <- ssa(modelTs, L = 25, column.projector='centering', row.projector='centering')
plot(modelTsSSA, type = "vectors")

plot(modelTsSSA, type = "paired")

modelTsSSASignal <- reconstruct(modelTsSSA, groups = list(1:2))

plot(modelTs, type = "l")
lines(modelTsSSASignal$F1, col = "red")

Toeplitz SSA

Теплицев SSA применяется как улучшение basic-SSA для стационарных рядов.

x <- 1:100
modelTs <- 5*cos(2*pi*(1/25)*x) + 2*cos(2*pi*(1/5)*x) + rnorm(50)

plot(modelTs, type="l")

modelTsSSA <- ssa(modelTs, kind = "toeplitz-ssa")

plot(modelTsSSA, type = "vectors")

plot(modelTsSSA, type = "paired")

modelTsSSASignal <- reconstruct(modelTsSSA, groups = list(1:4))

plot(modelTsSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))

plot(modelTs, type = "l")
lines(modelTsSSASignal$F1, col = "red")

Прогнозы: рекурентный и векторный

02.05

Рекурентный прогноз - rforecast

forcastR <- rforecast(tsSSA, groups = list(trend = c(1:4, 7:10), season = c(5, 6, 12, 13)), len = 12, only.new = F)

plot(unempl.male, type="l")
lines(data.frame(forcastR)$trend, col="red")
legend("topleft", legend = c("Real TS","RForcast trend"), col = c("black","red"), lty=1)

C отсечением периода

Применяем SSA к ряду с отсеченным периодом

tsCut = ts(unempl.male[1:(length(unempl.male)-12)], start=c(1948, 1), frequency=12)

tsCutSSA <- ssa(tsCut,  L = length(tsCut) %/% 2)
plot(wcor(tsCutSSA, groups = 1:30))

plot(wcor(tsCutSSA, groups = 1:15))

plot(tsCutSSA, type="vectors", idx = 1:20)

tsCutSSASignal = reconstruct(tsCutSSA, groups <- list(c(1:3, 6:10, 14:15), c(4,5,12,13)))

plot(tsCutSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))

forcastRCut <- rforecast(tsCutSSA, groups = list(trend = c(1:3, 6:10, 14:15), season = c(4,5,12,13)), len = 12, only.new = F)
plot(c(tsCut), type="l")
lines(c(data.frame(forcastRCut)$trend), col="red")
lines(c(data.frame(forcastR)$trend), col="blue")
legend("topleft", legend = c("Cut TS","RForcast cut trend", "RForcast full trend"), col = c("black","red", "blue"), lty=1)

Векторный прогноз - vforecast

forcastV <- vforecast(tsSSA, groups = list(trend = c(1:4, 7:10), season = c(5, 6, 12, 13)), len = 12, only.new = F)

plot(cbind(unempl.male, forcastV$trend), plot.type = "single", 
     col = c("black", "red"), ylab = NULL)

legend("topleft", legend = c("Real TS","VForcast full trend"), col = c("black","red"), lty=1)

C отсечением периода

forcastVCut <- vforecast(tsCutSSA, groups = list(trend = c(1:3, 6:10, 14:15), season = c(4,5,12,13)), len = 12, only.new = F)


plot(c(tsCut), type="l")
lines(c(data.frame(forcastVCut)$trend), col="red")
lines(c(data.frame(forcastV)$trend), col="blue")
legend("topleft", legend = c("Real TS","VForcast cut trend", "VForcast full trend"), col = c("black","red", "blue"), lty=1)

Восстановление ряда по ЛРФ

25.04/02.05

Ряд 1: Экспонента

library(Matrix)

Находим ранг матрицы для экспоненциального ряда:

N = 10
LRF = 3.1^(1:N)
rankMatrix(hankel(LRF))[1]
## [1] 1

Ожидаемо, он равен одному. Далее находим сигнальные корни ряда:

LRF.ssa = ssa(LRF, L=2, method = "svd")
l = lrr(LRF.ssa, groups = list(1))
plot(l)

roots_my = roots(l)
K = length(roots_my)
S_n = LRF[1:K]
lin_sys = data.frame(S_n)
vars = matrix(nrow = K, ncol = K)
for(i in 1:K){
  vars[i]=roots_my[i]^(1:K)
}
lm0 = lm(formula = S_n ~ 0 + ., data = data.frame(vars))
print(lm0)
## 
## Call:
## lm(formula = S_n ~ 0 + ., data = data.frame(vars))
## 
## Coefficients:
## vars  
##    1

Находим значения по линейной рекурентной формуле:

pred_S_n = lm0$coefficients[[1]]*(roots_my[1]^(1:N))

print(LRF)
##  [1]     3.1000     9.6100    29.7910    92.3521   286.2915   887.5037
##  [7]  2751.2614  8528.9104 26439.6222 81962.8287
S_n = LRF[1:N]
lin_sys = data.frame(S_n)
vars = matrix(nrow = N, ncol = K)
for(i in 1:K){
  vars[i]=roots_my[i]^(1:N)
}
lm0 = lm(formula = S_n ~ 0 + ., data = data.frame(vars))
print(lm0)
## 
## Call:
## lm(formula = S_n ~ 0 + ., data = data.frame(vars))
## 
## Coefficients:
## vars  
##    1

Сравниваем исходный ряд с полученным

pred_S_n = lm0$coefficients[[1]]*(roots_my[1]^(1:N))

print(LRF)
##  [1]     3.1000     9.6100    29.7910    92.3521   286.2915   887.5037
##  [7]  2751.2614  8528.9104 26439.6222 81962.8287
print(pred_S_n)
##  [1]     3.1000     9.6100    29.7910    92.3521   286.2915   887.5037
##  [7]  2751.2614  8528.9104 26439.6222 81962.8287

Ряд 2: Линейная функция

N = 10

LRF = 2*(1:N)+5

Находим ранг матрицы:

rankMatrix(hankel(LRF))[1]
## [1] 2

Линейная функция задается двумя компонентами.

LRF.ssa = ssa(LRF, L=3, method = "svd")
l <- lrr(LRF.ssa, groups = list(1:2))
roots_my <- roots(l)
print(l)
## [1] -1  2
## attr(,"class")
## [1] "lrr"
print(roots_my)
## [1] 1.0000002 0.9999998

Получаем два комплексно-сопряженных корня

print(2 * pi / Arg(roots_my))
## [1] Inf Inf
print(Mod(roots_my))
## [1] 1.0000002 0.9999998
roots_my2 <- rep(Re(mean(roots_my)), length(roots_my))
print(roots_my2)
## [1] 1 1

Находим коэффиценты исходного ряда

m = length(roots_my2)
s_n = LRF[1:N]
vars <- matrix(nrow = N, ncol = m)
for (i in 1:m) {
  vars[ ,i] <- (1:N) ^ (i-1) * roots_my2 ^ (1:N)
}

lm0 <- lm(s_n ~ 0 + ., data = data.frame(vars))
print(lm0)
## 
## Call:
## lm(formula = s_n ~ 0 + ., data = data.frame(vars))
## 
## Coefficients:
## X1  X2  
##  5   2

Сравниваем исходный ряд с полученным

pred_S_n = lm0$coefficients[[1]][1] * (vars[,1]) + lm0$coefficients[[2]][1] * (vars[,2])
print(LRF)
##  [1]  7  9 11 13 15 17 19 21 23 25
print(pred_S_n)
##  [1]  7  9 11 13 15 17 19 21 23 25

Восстановление ЛРФ реального временного ряда

Проанализируем собственные числа ряда

tsSSA$sigma
##  [1] 367075.857  74505.552  63304.591  35713.403  34728.609  34583.735
##  [7]  27721.976  22861.296  22230.219  18781.365  18634.256  18365.270
## [13]  18171.464  14245.491  13430.646  13113.963  10353.466  10239.538
## [19]   9456.080   9280.711   9009.749   8960.235   8234.105   8146.040
## [25]   7272.790   6845.598   6199.300   6024.433   5492.880   5475.516
## [31]   5071.906   4194.408   4169.282   3688.478   3611.481   3376.238
## [37]   3365.972   3203.616   3192.859   3179.853   3072.833   3014.386
## [43]   2996.338   2987.268   2984.135   2906.750   2672.857   2504.003
## [49]   2419.366   2386.851
plot(tsSSA$sigma, type="l")

Видим 13 “ненулевых” собственных числа.

rk = 13

par = parestimate(tsSSA, groups = list(1:rk), method = "esprit")

modulusReal = par$moduli
periodsReal = par$periods
o = order(abs(periodsReal), decreasing = TRUE)
r0 = reconstruct(tsSSA, groups = list(signal = 1:rk))$signal


len <- rk
vars <- matrix(nrow = len, ncol = rk)

for (i in 1:rk) {
  if (periodsReal[i] == Inf)
    vars[, i] <- modulusReal[i]^(1:len)
  else if (periodsReal[i] == 2)
    vars[, i] <- (-modulusReal[i])^(1:len)
  else if (periodsReal[i] > 0)
    vars[, i] <- 
      modulusReal[i]^(1:len) * sin(2 * pi * (1:len) / periodsReal[i])
  else
    vars[, i] <- 
      modulusReal[i]^(1:len) * cos(2 * pi * (1:len) / periodsReal[i])
}
lm0 <- lm(r0[1:len] ~ 0 + ., data = data.frame(vars))
coefs0 <- coef(lm0)
print(round(coefs0), digits = rk)
##        X1        X2        X3        X4        X5        X6        X7        X8 
##  -4965469  -5555248   7213618  54636156 -48963499      1651       739       205 
##        X9       X10       X11       X12       X13 
##       -65   1225627   -117262        NA        NA
idx <- seq(1, rk)
coefs.c.phase <- numeric(length(idx))
phases.c <- numeric(length(idx))
periods.c.phase <- numeric(length(idx))
moduli.c.phase <- numeric(length(idx))
for (i in seq_along(idx)) {
  
  periods.c.phase[i] <- periodsReal[idx[i]]
  moduli.c.phase[i] <- modulusReal[idx[i]]
  coefs.c.phase[i] <- coefs0[idx[i]]
  phases.c[i] <- atan2(coefs0[idx[i] + 1], coefs0[idx[i]])
  if(i == rk){
    phases.c[i] <- atan2(coefs0[idx[i]], coefs0[idx[i]])
  }
}

s = data.frame(periods = periods.c.phase, phases = phases.c, 
                 coefficients = coefs.c.phase, 
                 moduli =moduli.c.phase)
print(s)
##        periods      phases  coefficients    moduli
## 1    62.059089 -2.30019430 -4.965469e+06 1.0102164
## 2   -62.059089  2.22704139 -5.555248e+06 1.0102164
## 3   207.578092  1.43952547  7.213618e+06 1.0031497
## 4  -207.578092 -0.73069723  5.463616e+07 1.0031497
## 5          Inf  3.14155894 -4.896350e+07 1.0007054
## 6    11.995011  0.42101827  1.650952e+03 0.9986643
## 7   -11.995011  0.27033942  7.392869e+02 0.9986643
## 8     6.007629 -0.30780195  2.048739e+02 0.9982437
## 9    -6.007629  1.57084947 -6.513055e+01 0.9982437
## 10   38.009961 -0.09538464  1.225627e+06 0.9969873
## 11  -38.009961          NA -1.172619e+05 0.9969873
## 12   50.401313          NA            NA 0.9907013
## 13  -50.401313          NA            NA 0.9907013

ARIMA

Автоковариационная функция ACF

ts1 <- read.csv(file = "ts1.txt", header = TRUE, as.is = FALSE)
ts5 <- read.csv(file = "ts5.txt", header = TRUE, as.is = FALSE)

plot(ts1$x,  type = "l")

plot(ts5$x,  type = "l")

acf(ts1, main = "Функция для данных ts1")

acf(ts5, main = "Функция для данных ts5")

Функция PACF

По количеству столбцов, значимо отличающихся от нуля, можно определить количество параметров в авторегрессии.

pacf(ts1, main = "Функция для данных ts1")

pacf(ts5, main = "Функция для данных ts5")

Оценка параметров модели

fit1 <- auto.arima(ts1)
fit5 <- auto.arima(ts5)


fit1
## Series: ts1 
## ARIMA(2,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ar2     ma1     ma2     mean
##       -0.2652  0.6067  1.1344  0.1612  99.9179
## s.e.   0.0481  0.0335  0.0542  0.0454   0.1132
## 
## sigma^2 = 1.064:  log likelihood = -1448.07
## AIC=2908.15   AICc=2908.23   BIC=2937.59
fit5
## Series: ts5 
## ARIMA(3,0,0) with zero mean 
## 
## Coefficients:
##          ar1      ar2     ar3
##       0.3103  -0.1605  0.0984
## s.e.  0.0222   0.0230  0.0223
## 
## sigma^2 = 0.3961:  log likelihood = -1910.39
## AIC=3828.78   AICc=3828.8   BIC=3851.18
forecastedValues1 <- forecast(fit1, 20)
forecastedValues5 <- forecast(fit5, 40)

plot(forecastedValues1, main = "Graph with forecasting FIT1",
col.main = "darkgreen") 

plot(forecastedValues5, main = "Graph with forecasting FIT5",
col.main = "darkgreen") 

ARIMA на реальных данных + (ARIMA VS SSA VS ETS)

fitArima <- auto.arima(unempl.male)
checkresiduals(fitArima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,2)(2,1,1)[12] with drift
## Q* = 20.36, df = 17, p-value = 0.2562
## 
## Model df: 7.   Total lags used: 24
summary(fitArima)
## Series: unempl.male 
## ARIMA(2,0,2)(2,1,1)[12] with drift 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     sar1    sar2     sma1   drift
##       1.6149  -0.6340  -0.4319  0.1149  -0.0108  0.1072  -0.8538  5.8252
## s.e.  0.1395   0.1366   0.1483  0.0671   0.1022  0.0932   0.0909  3.6365
## 
## sigma^2 = 16160:  log likelihood = -2484.99
## AIC=4987.98   AICc=4988.44   BIC=5023.81
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.6993681 123.9684 91.59157 -0.34855 5.249927 0.1989636
##                      ACF1
## Training set 0.0005744349
fitArimaForecast <- forecast(fitArima, h = 24)
autoplot(fitArimaForecast)

tsSSABforecast <- bforecast(tsSSA, groups = list(c(1:4, 7:10)), len = 24, R = 50, only.new = F)


tsSSAETS <- hw(ts(unempl.male[1:(length(unempl.male)-24)], start=c(1948, 1), frequency=12), seasonal="multiplicative", h = 24)

autoplot(tsSSABforecast)

plot(tsSSAETS)

RMSE SSA vs ARIMA vs ETS

sd(unempl.male)
## [1] 794.9778
tsTest = tail(unempl.male, 24)

tsSSABforecast <- bforecast(tsSSA, groups = list(c(1:4, 7:10)), len = 24, R = 50, only.new = T)

mseSSA = mean((data.frame(tsSSABforecast)$Value - c(tsTest))**2)

mseARIMA = mean((c(fitArimaForecast$mean) - c(tsTest))**2)

mseETS = mean((c(tsSSAETS$mean) - c(tsTest))**2)

RMSE SSA

print(sqrt(mseSSA))
## [1] 1399.446

RMSE ARIMA

print(sqrt(mseARIMA))
## [1] 1024.97

RMSE ETS

print(sqrt(mseETS))
## [1] 889.6957

SSA и разладка

library(Rssa)
library("lattice")
firstHalf = c(10*sin(2*pi*(1:72)/6))
secondHalf = c(30*sin(2*pi*(73:144)/6))
tsdestruct = c(firstHalf,secondHalf)
plot(tsdestruct, type = "l", col = "red")

s <- ssa(tsdestruct, L = 12)
w <- wcor(s, groups = 1:10)
plot(w)

plot(s, type = "vectors", idx=1:10)

r <- reconstruct(s, groups = list(Trend = c(1)))

tsdestruct_res <- residuals(r)
N <- length(tsdestruct_res)
rank <- 2
periods <- function(M, L) {
  ts(sapply(1:(N - M),
            function (i) {
              s <- ssa(tsdestruct_res[i:(i + M - 1)], L = L)
              par <- parestimate(s, groups = list(c(1:rank)), 
                                 method = "esprit")
              abs(par$periods[1])
            }),
     start = time(tsdestruct)[M + 1], delta = 1)
}
per12 <- periods(12, 6)
per24 <- periods(24, 12)
lattice::xyplot(plot.method = "xyplot",per12, type = "l")

M <- 12; L <- M / 2
hm <- hmatr(tsdestruct_res, B = M, T = M, L = L, neig = rank)
plot(hm)